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ABSTRACT 

We present results from a numerical study of the runaway instability of thick discs 
around black holes. This instability is an important issue for most models of cosmic 
gamma-ray bursts, where the central engine responsible for the initial energy release 
is such a system consisting of a thick disc surrounding a black hole. We have carried 
out a comprehensive number of time-dependent simulations aimed at exploring the 
appearance of the instability. Our study has been performed using a fully relativistic 
hydrodynamics code. The general relativistic hydrodynamic equations are formulated 
as a hyperbolic flux-conservative system and solved using a suitable Godunov-type 
scheme. We build a series of constant angular momentum discs around a Schwarzschild 
black hole. Furthermore, the self-gravity of the disc is neglected and the evolution of 
the central black hole is assumed to be that of a sequence of exact Schwarzschild black 
holes of varying mass. The black hole mass increase is thus determined by the mass 
accretion rate across the event horizon. In agreement with previous studies based on 
stationary models, we find that by allowing the mass of the black hole to grow the 
disc becomes unstable. Our hydrodynamical simulations show that for all disc-to-hole 
mass ratios considered (between 1 and 0.05), the runaway instability appears very fast 
on a dynamical timescale of a few orbital periods, typically a few 10 ms and never 
exceeding 1 s for our particular choice of the mass of the black hole (2.5 M©) and a 
large range of mass fluxes (m>10 -3 M Q /s). The implications of our results in the 
context of gamma-ray bursts are briefly discussed. 

Key words: accretion, accretion discs - black hole physics - hydrodynamics - insta- 
bilities - gamma rays: bursts. 



1 INTRODUCTION 

Thick accretion discs are probably present in many astro- 
physical objects, e.g. quasars and other active galactic nu- 
clei, some X-ray binaries and microquasars, and the cen- 
tral engine of gamma-ray bursts (GRBs hereafter). They 
have been s tudied in great detail by many authors (see e.g. 
Rees (1984 ) and references therein). In particular, it is well 
known that in a system formed by a black hole surrounded 
by a thick disc, the gas flows in an effective (gravitational 
plus centrifugal) potential, whose structure is comparable 
with that of a close binary. The Roche torus surrounding 
the black hole has a cusp-like inner edge located at the La- 
grange point Li where mass transfer driven by th e radial 
pressure gradient is possible (Kozlowski et al. 1978). 



The so-called runaway instability of suc h systems was 
first discovered by Abramowicz et al. (1983 ). The underly- 



ing mechanism is the following: due to accretion of mate- 
rial from the disc, the mass of the black hole increases and 
the gravitational field of the system changes. Therefore, an 
accretion disc can never reach a completely steady state. 
Starting from an initial disc filling its Roche lobe so that 
mass transfer is possible through the cusp located at the Li 
Lagrange point, two evolutions are feasible when the mass 
of the black hole increases: (i) either the cusp moves inwards 
towards the black hole, which slows down the mass transfer, 
resulting in a stable situation, or (ii) the cusp moves deeper 
inside the disc material. In this case the mass transfer speeds 
up, leading to the runaway instability. 



In their first study, Abramowicz et al. (1983) analyzed 



the effect of the mass transfer under many simplifying as- 
s umptions (a pseudo-New tonian potential for the black hole 
(Paczyhski & Wiita 1980), constant angular momentum in 



the disc, and approximate treatment of the disc self-gravity) . 
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Table 1. Summary of representative results concerning the runaway instability on stationary models of thick discs around black holes. 



Reference 



Framework 



Angular momentum BH rotation Disc self-gravity Results 



Abramowicz et al. 



Wilson (198W 



Khanna fc Chakrabarti (1992) 



Nishida ct al. (1996 ) 



Daigne fc Moclikovitch ( 1997) 



Abramowicz et al. (1 996 ) 



Masuda ct al. ( 1998)' 



Lu et al. (2000) 



Pseudo-Newt. 
GR 

Pseudo-Newt. 
GR 

Pseudo-Newt. 
GR 

Pseudo-Newt. 
Pseudo-Newt. 



= constant no 

= constant yes 

= constant yes 

= constant yes 

oc ro Q , < a < 0.5 no 

oc ra Q , < a < 0.5 yes 

oc zu a , a = 0.2 no 

oc ro a , < a < 0.5 no 



approximate 

neglected 

approximate 

exact 

neglected 

neglected 

exact 

neglected 



unstablc 1 

stable 

unstable 

unstable 2 

stable 3 

stable 4 

unstable 5 

stable 6 



Notes : 

1 for a large range of disc-to-hole mass ratio and disc inner radius. 

2 study made for only four intial models. 

3 The mass of the bla ck hole is 2.44 Mm. The mass of the disc is 0.36 M c . 

4 Qamo naramoti: 



Same parameters as 



The system is stable for a > o cr with a CI ~ 0.1. 
r decreases when the black hole rotation increases. 
Mochkovitch [1997 ). The system becomes unstable for large transfer of mass only. 



Daigne < 


k Mochkovitch ( 1997 


Daigne < 


fe Mochkovitch (1997 


Daigne & 


- Mochkovitch (1997|) 



Same conclusion as Daigne fc Mochkovitch (1997f for a completely different range of masses (massive black hole of mass 10 6 Mq) 



These authors found that the runaway instability occurs for 
a large range of parameters such as the disc-to-hole mass ra- 
tio and the location of the disc inner radius. More detailed 
studies followed, whose main conclusions are summarized 
in Table [j]. Notice that all these studies are based on sta- 
tionary models in which a fraction of the mass and angular 
momentum of an initial disc filling its Roche lobe is trans- 
fered to the black hole, and the new gravitational field is 



used tq compute the new position ot the cusp, which con- 
trols whether the runaway instability occurs or iiuL The 

conclusions of these studies have yet to be confirmed on a 
dynamical framework. 



From Table hi one sees that (i) taking intc 



self-grakjty oXlhg disc sggmg ta foggr the instability f|Khanna 




aunt the 



& Chajcrabarti 1992|; |Nishida et al. 1996 
1996 
effect 



Masuda et al 



the rotat i on of the black hole has a stabilizing 
ilson 1984; Abramowicz et al. 1998); (iii) taking 



into account a non-constant distribution of the angular mo 
ment um (increasing outwards) has a strong stabilizing e f- 



fect (Daigne fc Mochkovitch 1997; Abramowicz et al. 1998); 



(iv) using a fully relativistic potential instead of a pseudo- 
Newtonian potential for the black h ole seems to act in the 
direction of favoring the instability (Nishida et al. 1996). It 
also becomes evident from Table hi that there is not still a 
final consensus about the very existence of the ins tability. 
In the fully relativistic study of Nishida et al. (1996), it was 



shown that the runaway instability occurs when the angu- 
lar momentum is constant in the disc. However, the work of 



Daigne & Mochkovitch (1997) showed the stabilizing effect 



fer is supposed to diverge, so that the thick disc could, in 
principle, be completely destroyed. 



In particular, Nishida et al. (1996) pointed out that the 
runaway instability could be a severe problem for most cur- 
rent models of GRBs. These models usually assume that the 
central engine is such a system consisting of a black hole and 
a thick disc, either formed at the late stages of t he coales 



cence of two neutron stars ( Kluzniak fc Lee 1998 



Ruffert & 



Janka 199!: ; ghibata fc Uryti 2000|) or after t he gravitationa l 



collapse of a massive star (Paczynski 1986; Woosley 1993) 
Notice that in the latter case, the situation is somewhat 
more complicated because the mass of the disc increases 
with time, as l ong as the collap se proceeds ( MacFadyen & 
Woosley 1999| ; [Aloy et al. 200C| ). The energy which can be 



extracted from this system comes from two reservoirs: the 
energy released by the accretion of disc material on to the 
black hole and the rotational energy of the black hole itself, 
which can be extracted via the Blandford-Znajek mecha- 
nism (Blandford fc Znajck 1977). This amount of energy 



(a few 10 53 -10 54 erg depending on the disc mass and the 
black hole rotation and mass) is sufficient to power a GRB 
if the energy released can be eventually converted into j- 
rays with a large efficiency of about a few percent. This 
number, however, depends strongly on the central engine 
model and on the expected beaming of the outflow, which is 
currently very poorly constrained. Such conversion cannot 
be done close to the source, the luminosity considered here 
being many orders of magnitude larger than the Edding- 
ton luminosity of the system. The energy is first injected 
into a very optically thick wind. This wind is accelerated 
via some mechanism which is still unknown but which in- 



of a distribution of angular momentum in the disc increasing 
outwards. It is worth pointing out that the complete calcu- 

lation in this setup, i.o,, in general relativity and including volves probably MHD processes ( [Thompson 1994 



Self-graidty and a rotating black hnlp is pvtrpmely rnmpW 

and has not been done yet. 

The consequences of the runaway instability could be 
very important in many cases. First because this is a purely 
dynamical effect, so that its time-scale is extremely short 
(only a few ms for a stellar mass black hole) . This means 



Meszaros 



Rees 1997b) and becomes eventually relativistic. The ex- 



istence of such a relativistic wind has been directly inferred 
from the observations of radio scintillation in GRB 970508 



(Waxman et al. 1998) and is also needed to solve the so- 



that this instability should happen before any other pro- 



called compactness problem and avoid photon-photon an- 
nihilation along the line of si ght. Averaged Lorentz factors 
larger than 100 ar e required ( Baring fc Harding 1997 ; Lith- 



wick 



cesses (like the viscous transport of angular momentum) 
could play a role. Secondly, because the radial mass trans- 



Sari 2001 ). The observed emission is produced at 
large distances from the source (r > 10 11 -10 12 cm), proba- 
bly via the formation of shock waves, either within the wind 
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itself (inte rnal shock model, proposed for the prompt 7-ray 
emission, (|Rccs fe Mcszaros 1994; Kobayashi ct al. 19 97^ 



numerical methods for their solution (see e.g. Font (200C ) 



Daigne|fc Mochkovitch 1998[)) or caused by the deceleration 
of the wind by the external medium (external sho ck model 



which r eproduces correctly the afterglow emission ( Meszaroi 
& Rees 1997a; 3ari et al. 



199J)) 



The above general scenario clearly presupposes that the 
black hole plus thick disc system is stable enough to survive 
for a few seconds (in particular, the internal shock model 
implies that the duration of the energy release by the source 
has a duration comparable with the observed duration of 
the GRB). However, if the runaway instability would occur, 
the disc would fall into the hole in just a few milliseconds! 

In order to establish the nature of instabilities in ac- 
cretion flows one must rely upon highly accurate, time- 
dependent, nonlinear numerical simulations in black hole 
spacetimes. These simulations are scarce and they have been 
mostly performed using a Newtonian or a pseudo-Newtonian 
poten tial which mimics the exist ence of an "event hori- 
on" riPaczvnski Wiit.a. 1 98(t fsee iFWum of al. f1988h:lpTl 
paloizo 1 fc Szuszkiewicz (19941); h^uinenslichfiv_e^ > aL_^199£ \ 



for simulations in the context of thick discs). Concerning 
relativistic simulations Wilson (1972) was the first to study 
numerically the time-dependent accretion of matter on to 
a rotating black hole. His simulations showed the forma- 
tion of thick accretion discs. In a subsequent work Hawley 
et al. ( |l984a| Jb|) studied the evolution and development of 
nonlinear instabilities in pressure-supported accretion discs 
formed as a consequence of the spiraling infall of fluid with 
some amount of angular momentum. Their constant angu- 
lar momentum initial models were computed fol lowing the 
analyti c the ory of relativistic discs developed by Kozlowski 
et al. (l|978|). The code developed by |Hawley et al. (1984ap |) 
was capable of keeping stable discs in equilibrium as well as 
of following the fate of initially unstable models. Yokosawa 
(1995^tudied the structure and dynamics of relativistic ac- 
cretion discs and the transport of energy and angular mo- 
mentum in magnetohydrodyn amical accretion on to a rotat 



ing black hole. More recently, [gumenshchev fc Beloborodov 
(1997^>erformed a similar study as jiawtey*^t*!iL*7T984a| ,}3[) 



including the rotating black hole case and using improved 
numerical methods based on Riemann solvers. They found 
that the structure of the innermost part of the disc depends 
strongly on the black hole spin and, at the same time, they 
were able to confirm numerically the expected analytic de- 
pendence of the mass accretion rate with the gravitational 
energy gap at the cusp of the torus. 

A time-dependent and fully relativistic study of the run- 
away stability has not ye t been presented in the literature. 



Masuda & Eriguchi (1997) performed a time-dependent sim- 



ulation in a pseudo-Newtonian framework, using an SPH 
method so that the self-gravity was taken into account. Our 
work aims at providing the first relativistic study. We have 
investigated the likelihood of the instability in thick discs 
of constant angular momentum. In this first investigation 
we present results for a Schwarzschild black hole. The rotat- 
ing case and the inclusion of accretion discs with a distri- 
bution of angular momentum increasing outwards will be 
presented elsewhere. Our simulations are performed in a 
fully relativistic framework, using a 3+1 conservative for- 
mulation ofoth^Jrvdxod^mamic equations on a curved back- 
ground ( Banyuls et al. 1997 ) and employing Godunov-type 



and re ferences therein). As in the work of Hawley et al 
( 1984a ,^) and [gumenshchev fc Beloborodov (1997| ) we ne- 



glect viscous and radiative processes assuming that the flow 
is isentropic. This is justified as we are only interested in phe- 
nomena occurring on a dynamical timescale. Furthermore, 
as a major simplification of the computational burden, we 
neglect the self-gravity of the disc. The spacetime dynam- 
ics has hence been treated in a simple way, assuming that 
the background spacetime metric is nothing but a sequence 
of stationary exact black hole solutions (Schwarzschild or 
Kerr) of the Einstein equations, whose dynamics is governed 
by the increase of the black hole mass (and angular momen- 
tum). Such growth is controlled only by the rate at which 
the mass accretes on to the hole. We are aware that a fully 
consistent approach to study the problem would imply solv- 
ing the coupled system of Einstein and hydrodynamic equa- 
tions in three-dimensions. This seems still a daunting task 
despite some major recent progress on numeric ally stable 
integrations of the Einstein eq uations (see, e.g. 



et al. (2000); Font et al. (2001) and references therein). We 



Alcubierre 



believe that, as a first step, our approach is justified and can 
provide some insight on the problem. 

The paper is organized as follows: Section || reviews 
briefly the theory of stationary, relativistic accretion discs of 
constant angular momentum. In Section ^ we introduce the 
equations of general relativistic hydrodynamics in the way 
they are used in our code. The numerical schemes we use 
to solve those equations, as well as other relevant aspects 
of our numerical code are described in Section ^. In Sec- 
tion ^| we test the capability of the code in keeping numeri- 
cally the equilibrium of stationary models in time-dependent 
simulations. The main results of the paper are presented in 
Section [] which contains the simulations of the runaway 
instability. Finally, Section [j] presents a summary of our in- 
vestigation. 

Unless explicitly stated we are using geometrized units 
(G — c — 1) throughout the paper. Greek (Latin) indices run 
from to 3 (1 to 3). We also use the signature -+++. Usual 
cgs units are obtained by using the gravitational radius of 
the black hole, r g = GMbh/c 2 , as unit of length. Although 
this paper is only concerned with the Schwarzschild case we 
present all expressions for the most general setup of a Kerr 
black hole in order to refer to them in the forthcoming papers 
of this series. The code is already written for the Kerr metric 
and the results here presented have been obtained by just 
setting to zero the black hole angular momentum parameter. 



2 STATIONARY MODELS OF THICK 
ACCRETION DISCS 

The theory of stationary, relativistic thick discs ( or tori) of 
constant angular m omentum was first derived by Fishbone 
&i Moncrief (1976) for isentropic discs and by 



Kozlowski 

et*lu"**(T978| ) for discs obeying a barotropic equation of state 
(EoS). Although this theory is by now well-established and 
has been presented in great detail elsewhere, we have chosen 
to include here the relevant aspects and definitions in order 
to facilitate the self-consistency of the paper. 
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2.1 Definitions 

2.1.1 Assumptions 

In order to build up constant angular momentum configura- 
tions for a given metric we assume a number of conditions. 
Firstly, the torus is stationary and axially symmetric. We 
adopt the standard spherical coordinates (t, r, 9, (f>) in which 
the metric coefficients neither depend on the time coordinate 
t (stationarity) or on the azimuthal coordinate <j> (axisym- 
metry). The line element of the spacetime is given by 

ds 2 = gudt 2 + 2guj,dtd<t> + g^dcp 2 + g rr dr 2 + geedQ 2 . (1) 



We adopt the following convention: the angular momentum 
of the black hole is positive and the matter of the disc ro- 
tates in the positive (negative) direction of <f> for a prograde 
(retrograde) disc. We define the following quantity which is 
a relativistic generalization of the "distance to the axis" : 



: 9t<t> — gag,/,,/,, 



(2) 



which in the Newtonian limit is simply zu = rsmQ. 

The EoS of the fluid is assumed to be barotropic, so 
that p = p(e) where p is the pressure and e is the total 
energy density. We define p as the rest mass density, w = 
e + P as the enthalpy and h — w/ p as the specific enthalpy. 
The four velocity of the fluid is it M = (u*,O,O,w ) with the 
normalization condition u 2 = u u u^ = — 1. 



2.1.2 Von Zeipel's cylinders 

The Lagrangian angular momentum (angular momentum 
per unit inertial mass) and the angular velocity are respec- 
tively defined by 

1 = -?* 



ut 



and 



u 

From u u 



g^v^ it follows that 



9t4> + 9<t>(t>Q 



and Q. - 



gt</> + gttl 



+ 



(3) 



(4) 



(•») 



gtt + gt<p£l 

For a barotropic EoS, the equi-Z and equi-fi surfaces coin- 
cide and they are called the von Zeipel's cyl inders. Their 
cylind er-like topology has been proved by Abramowicz 
(1974). If the metric is known, Eq. allows to construct 



the von Zeipel's cylinder defined by Iq and fio by solving the 
following equation 



gttlo + 9t<t> (1 + fio'o) + g$<t>Qo = 



(6) 



In particular, if the distribution of the angular momentum 
I = lcq(r) is given in the equatorial plane 9 — n/2, then 
the corresponding distribution of the angular velocity in the 
equatorial plane is given by 



g t 4,(r,n/2)+gtt{r,n/2)l(r) 
g<M.(r, 7r/2) + g t <t,{r, n/2)l(r) 



(7) 



The equation of the von Zeipel's cylinder intersecting the 
equatorial plane at a given radial point r — ro is given by 
Eq. (|) so that 

l 2 {ro) [gtt(r, 9)g t 4,(r , tt/2) - g t<t ,(r, 9)g u (r , n/2)] 



+l{ro) [gtt(r, 9)g 4>4 ,(r , tt/2) - g^(r, 9)g tt (r , tt/2)] 

+ {gtj>(r, 9)g tj ,j,(r , tt/2) ~ g^{r,9)g t $(r , tt/2)] = 0. 

(8) 

Notice that in the case of the Schwarzschild metric, g t ^ = 0, 
and this equation reduces to 



gtt(r, 9)g^(r ,ir/2) - g^(r, 6)g u (ra, tt/2) = 



(9) 



which is independent of the distribution of angular momen- 
tum. 



2.1.3 Equipotentials 

The dynamics of the gas flow is governed by the relativistic 
Euler equation, whose integral form is 

rl 



W-W in = In (-«*) - In (- u t in ) - 



Qdl 



(10) 



The subscript "in" refers to the inner edge (in the equatorial 
plane) of the disc, where the pressure vanishes. The quantity 
W is defined by 



M -H ; „ = - f * 

w 



(11) 



In the Newtonian limit, the quantity W is the total (cen- 
trifugal plus gravitational) potential and Eq. (|ll|) is the in- 
tegral form of the equation of hydrostatic equilibrium. If the 
spacetime metric is known and if the distribution of angular 
momentum in the equatorial plane is given, so that the von 
Zeipel's cylinders have been computed, providing us with 
the value of I and Q at any given point within the disc, then 
the equipotentials surfaces W(r, 9) are easily computed from 
Eq. (|l(j) taking into account that ut can be expressed (from 
the 4- velocity normalization condition u 2 = — 1) as a func- 
tion of I by 



Ut 



g tt l 2 + 2g t<t> l + g^ 



(12) 



2.1.4 Construction of a thick disc 

In order to build a system consisting of a black hole sur- 
rounded by a thick disc we need several parameters: the 
mass M and the specific angular momentum a = J/M of 
the black hole, the distribution of the angular momentum 
of the disc in the equatorial plane Z eq (i~), the inner radius 
of the disc r ln and the EoS of the fluid material of the disc. 
The procedure is then the following: 

(i) Compute the metric coefficients in a Kerr background 
with free parameters a and M (these coefficients are given 
in the next Section). 

(ii) From the distribution of the angular momentum in 
the equatorial plane, solve Eq. (^|) to have the distribution 
of angular momentum l(r,9). Then the corresponding dis- 
tribution u t (r,9) can be evaluated from Eq. (p^[). 

(iii) From the value of the inner radius n n , compute the 
corresponding value of the angular momentum l ln = Z cq (ri n ). 
Then the function 



F(r,9) 



Qdl 

1 - m 



(13) 
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Figure 1. The von Zeipel's cylinders for the Schwarzschild met- 
ric. The critical cylinder (thick line) has a cusp located at r = 3 
in the equatorial plane. 



can be evaluated. 

(iv) Compute the potential 

W(r, 9) - W in = In (-u t (r, 9)) - In (- u t in ) - F{r, 9) . 

The constant Wi n is fixed by the convention that 
\im r ^ +ca W(r,9) = 0. 

(v) Compute all hydrodynamical quantities (p, p, e, w, 
etc) from the EoS and Eq. ([H]). In the following , we will 
only consider the particular case of isentropic fluids where 



v dp , h 

— = In- 
ly hi n 



(14) 



If the self-gravity of the disc is neglected, the procedure stops 
here. Otherwise, from the distribution of matter-energy we 
have computed, we need to solve the Einstein field equations 
to evaluate the new metric coefficients g^ v . Then the proce- 
dure starts again at step 2. Such cycles are repeated until 
convergence. 

Once a disc has been built, one can estimate its mass 
with the expression 



'Ti + T r r+ T! 



2?) sT~g dv 



(15) 



where T^ v is the stress-energy tensor which, in the case of a 
perfect fluid is defined by T M " = phu^u" + pg^ . In all the 
models we will consider in the following, we have p < e so 
that the mass of the disc can be accurately approximated 
by the rest-mass 



p [u<p u 
1-Q.l 



* - u t u) yf=g~ dV 
p^—g dV . 



(16) 



2.2 Constant angular momentum thick discs 
orbiting a Schwarzschild black hole 

In this subsection we describe in detail the particular case we 
have considered for the time-dependent simulations of this 
paper: the black hole is non-rotating (a Schwarzschild black 
hole), the self-gravity of the disc is neglected, the angular 
momentum / is constant and the disc obeys a polytropic 
EoS with 



P = Kp' 



(17) 



with k being the polytropic constant and 7 the adiabatic 
exponent. For commodity, we further assume that the mass 
of the black hole is equal to unity. The metric coefficients 
are given by 



9tt 

9<t><t> 



= -I1-- 

r 



(18) 

(19) 
(20) 

and the distance to the axis is simply zu 2 = r(r — 2) sin 2 9. As 
the black hole is non-rotating, we also consider only positive 
values of the angular momentum (prograde discs). 

The von Zeipel's cylinders are independent of I. Using 
Eq. the cylinder intersecting the equatorial plane at r = 
ro is given by 



(r - 2)r 3 



■rg(r-2) =0 



(21) 



The result is shown in Fig. EL Notice the cusp located at 
r = 3 in the equatorial plane (thick line). 

Since the angular momentum / is constant, the function 
F(r,9) vanishes. Using Eq. (|l2]), the component ut of the 4- 
velocity is given by 



Ut 



r sin 9 



(r - 2)P 



and the potential W reads 



I r s s 



: (r-2)sin 2 6i 
n 2 9-P{r- 2)' 



(22) 



(23) 



where we have used the condition W(r, 9) — + for r — ► +00 
to eliminate n n . 

We describe now the equipotentials which are either 
closed (W < 0) or open (W > 0). The marginal case W = 
is closed at infinity. The geometry of the equipotentials is 
fixed by the value of I. We impose that W(r) is defined 
everywhere outside the horizon in the equatorial plane, i.e. 
that the term r 3 — l 2 (r — 2) never vanishes. Then I < Z max = 
3^3 ~ 5.20. 

It can be easily shown that the presence of a cusp in 
the equatorial plane is related to the solutions of 



l K (r)=l 



(24) 



where ln{r) is the Keplerian angular momentum of a particle 
located at a radius r in the equatorial plane. It is given by 



lx(r) = 



(25) 



Two useful radii can be defined, the radius of the last 
(marginally) stable orbit r ms = 6 (corresponding to the 
minimum of lK(r)), and the radius of the last (marginally) 
bound orbit r m b = 4. The corresponding values of Ik are 
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Table 2. The possible equipotential configurations for a Schwarzschild black hole and a constant angular momentum disc. 



Angular momentum 


^"cusp 




W cusp 


^center 




^center 


Comments 


I < ims — 3.67 














No cusp. No center. Disc infinite. 


I = ims 2± 3.67 


^cusp — f*ms 


= 6 


< 


Reenter 


= r ms = 6 


< 


Cusp = center. Disc infinite. 


Us <C / < ^mb 


r mh < r cusp 




< 


Reenter 


^ fras 


< 


Cusp. Center. Disc closed. 


I = 'mb = 4 


'"cusp — r mh 


= 4 





Reenter 


~ 10.47 


< 


Cusp. Center. Disc closed at infinity. 


^mb *^ ^ ^ ^max 


^cusp < r mh 




> 


^'center 


> 10.47 


< 


Cusp. Center. Disc infinite* . 


' — ^max — 5.20 


^cusp = 3 




+oo 


Reenter 


~ 22.39 


< 


Cusp marginally defined. Center. Disc infinite* . 



Some closed equipotentials are still present around the center. 



Ims = ^ ~ 3.67 and l mb = 4. The properties of W(r, 9) 
are given in Table ^ for the different possible values of I. 
The corresponding equipotentials are drawn in Fig. ^ (right 
panel). The angular momentum and potential in the equato- 
rial plane are also drawn on the left panel of Fig. ^| It is clear 
that the most interesting case is case (3) where Z ms < I < lmb 
so that there exist a cusp, a center, and the equipotential 
of the cusp is closed. From Eq. ( p"7| ) and Eq. (jlij) we have 
indeed 



(26) 



W - Win = - In h = - In [ 1 + — ^— np 1 - 1 
V 7-1 

so that the matter can fill only the part where W < Wi n . The 
density and the pressure are then easily determined from h: 
, _ w in -w 



7 — 1 e 



w in -w 



7 

7 — 1 e 



K 

w ia -w 



i \ i— i 



1 



K 



(27) 
(28) 

(29) 



It is then possible to adjust the value of I and n n to fix the 
mass of the disc, which is given by 



, . r 3 sin 2 fl + (r-2)l 2 a . 

m ~ 2-7T lip = ; ; — r smOdOdr 

' ' V sin 2 - (r - 2)Z 2 



(30) 



If one iilnposes the condition that the disc is exactly filling its 



Roche lobe, the inner radius ri n = r cusp is fixed. Otherwise, 
instead of specifying r- m , one could prefer another parameter, 
such as the potential barrier (energy gap) at the inner edge 
defined as 



AW in = W in 



W a 



(31) 



The case AWi n < corresponds to a disc inside its Roche 
lobe. No mass transfer is possible. The case A Win > corre- 
sponds to a disc overflowing its Roche lobe: mass transfer is 
possible at the cusp. An analytic estimation for the mass flux 
(flux of rest mass density) was derived by Kozlowski et al 
(197&f}or this last case, showing the following dependence 



m oc AW2 



(32) 



3 HYDRODYNAMIC EQUATIONS IN A KERR 
BACKGROUND 

Our purpose is to evolve in time the initial data described 
in the previous Section. In order to do so we present in this 



Section the formulation of the general relativistic hydrody- 
namic equations in the form they have been implemented in 
our numerical code. 

Using Boyer-Lindquist (t, r, 9, (f>) coordinates, the Kerr 
line element, ds 2 = g^dx^dx" , reads 



ds 



-dt - 2a- 



i 2 . 2 , a 2 . 2 • 

+ -r-dr + g d& -\ -sm 



with the definitions: 



2Mr sin 2 9 

2 /ijj.2 



-dtd<f> 



A 

2 

Q 
S 



r 2 -2Mr + a 2 , 
r 2 + a 2 cos 2 6, 



(r + a ) — a A sin 9, 



(33) 



(34) 
(35) 
(36) 



where M is the mass of the black hole and a is the black hole 
angular momentum per unit mass (J/M). Notice that the 
geometrical factor g has not to be confused with the rest- 
mass density of the fluid, p. The above metric, Eq. (^), de- 
scribes the spacetime exterior to a rotating and non-charged 
black hole. The metric has a coordinate singularity at the 
roots of the equation A = 0, which correspond to the hori- 
zons of a rotating black hole, r = r± = M ± (M 2 — a 2 ) 1 ^ 2 . 
The "distance to the rotation axis" introduced in the previ- 
ous Section is given by zu 2 — g 2 ^ — gu9tj><j> = A sin 2 9. 

The {3+1} decomposition (see, e.g., Misner et al. 

(1973| )) of this form of the metric leads to a spatial 3- metric 



jij with non-zero elements given by 7 r 



Q /A, 790 



Q ,"/<t>4> = E/g sm 6. In addition, the azimuthal shift vector 



= gtcf, is given by 
2aMrsin 2 <9 

p 2 

and the lapse function is given by 
E 



(37) 



(38) 



The equations of general relativistic hydrodynamics are 
obtained from the local conservation laws of density current 
J* and stress-energy T MI/ 



with 

r = 



pan , u + pg , 



(39) 
(40) 

(41) 
(42) 
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Figure 2. Constant angular momentum thick disc orbiting a Schwarzschild black hole: Equipotcntials for different values of / (right 
column). The thick line corresponds to the equipotcntial of the cusp and a thick dot marks the center. The physically interesting case 
(corresponding to a thick torus) is case (3). 
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Figure 2 — continued Constant angular momentum thick disc orbiting a Schwarzschild black hole: Equipotentials for different values of 
I (right column). The thick line corresponds to the cquipotcntial of the cusp and a thick dot marks the center. 
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for a general EoS of the form p = p(p, e), e being the specific 
internal energy. The specific enthalpy is defined as h — 1 +e+ 
P/p- Furthermore, is the covariant derivative associated 
with the four- dimensional metric and u 11 is the fluid 4- 
velocity. The above expression of the stress-energy tensor 
corresponds to that of a perfect fluid. 



1 = 1 g t <t,,r + 1 gt<t>,e + 1 £ 



Following the general approach laid out in Banyuls et al 
(f 997^7|after the choice of an appropriate vector of conserved 
quantities, the general relativistic hydrodynamic equations 
can be written as a first-order flux-conservative hyperbolic 
system. In axisymmetry (d,/, = 0) and with respect to the 
Kerr metric such system adopts the form 



aU(w) + d(aF r (w)) + <9(aF a (w)) 



S(w) 



(43) 



dt dr 89 

In this equation the vector of (physical) primitive variables 
is defined as 



(44) 



where Vi(i = r,9,<f>) is the fluid 3- velocity, defined as v % = 
u z /au 1 + j3 % /a, with Vi = jijV 3 . On the other hand, the state 
vector (evolved quantities) in Eq. (143) is 



U(w) = (D, S r , Se,S^, t). 



(45) 



The explicit relations between the two sets of variables, U 
and w, are 

D = P T, 

S 3 = P hV 2 v J C7=r,M)> (46) 
r = phT 2 -p-D, 

with r being the Lorentz factor, Y = au* = (1 — n 2 ) _1//2 , 
with v 2 = jijv'vi . -phe S p ec i nc form of the fluxes, F l , and 
the source terms, S, read 



F r (w) = (Dv r , S r v r + p, Sgv r , S^v r , (r + p)v r ) , 

F"(W) = (Dv e ,SrV e ,SgV e +p,S tt> V <> , (T+ P )v e ) , 

(Si, 5*2, S3, S4, S5), 



S(w) 
with 

Si : 

s 2 ■ 

s 4 ■ 
s 5 ■ 



—ADv r — BDv 9 , 

-A(S r v r +p)- BSrv" + aC, 

-AS e v r - B(S e v e + p) + aD, 

—AS^v r — BScfyV + a£, 

—A(t + p)v r — B(t + p)v° + ajF, 



with the definitions 
A = 



r 1 E, r A 



EA , 



B = 
C = 

+ 

V = 



+ 



AE 



a I cot 9 sin 6 cos ( 

Q 2 



a 2 A 
2E 



sin 26» 



1 rrtrB 
;r + 1 g r 



g r r{T n Y r H + T rr T r rr + T ttH r r eg + T 4 " t T r 

2r^rL + 2T r9 r; e ), 



£0 

Tr6 . rrv 

aee, r + ± 



gee(T Vu + T T rr + T T ge +T vv Y (j>4 , 
2T" /> rL + 2T re F e rB ), 



(47) 
(48) 
(49) 

(50) 

(51) 
(52) 
(53) 
(54) 

(55) 
(56) 

(57) 
(58) 



- 2g t4> (T tT V\ r + T tB V\ e + r*Ti4 + T e *r^) 

- 2g 4 , 4> (T tr Ti + T t0 Ti + T^T% + T^Ti 4> ), (59) 

T = T tr a, r + T te a, 9 

- 2a(T tr V\ r + T^Tle + T r< *T^ + T^T^). (60) 

The ',' in the above expressions denotes partial differentia- 
tion and the Tt v stand for the Christoffel symbols (only the 
non-vanishing ones are displayed) which are obtained from 
the metric according to the usual definition 



(61) 



4 NUMERICAL METHOD 
4.1 The hydrodynamics code 

The hydrodynamics code used in our computations was orig- 
inally developed fo r studies of relativistic win d accretion 
on to black holes (Font & Ibaiiez (1998a b|) 




Font et al 



This code performs the numerical integration 
of system ([13[) using a so-called Godunov-type scheme. Such 
schemes are specifically designed to solve nonlinear hyper- 
bolic systems of conservation laws (see, e.g. Toro 1999 for 
definitions). In a Godunov-type method the knowledge of 
the characteristic structure of the equations is essential to 
design a solution procedure based upon either exact or ap- 
proximate Riemann solvers. These solvers compute, at ev- 
ery cell-interface of the numerical grid, the solution of local 
Riemann problems (i.e., the simplest initial value problems 
with discontinuous initial data). Therefore, they automat- 
ically guarantee the proper capturing of all discontinuities 
which may arise naturally in the solution space of a nonlin- 
ear hyperbolic system. 

The time update of system ( ji^ ) from t" to t n+1 is per- 
formed according to the following conservative algorithm: 



U" 



F'j-i 



/2,j) 



At Si j . 



(62) 



Index n represents the time level and the time (space) dis- 
cretization interval is indicated by At (Ar, A9). The numer- 
ical fluxes in the above equation, F r , F e are computed by 
means of the HLLE Riemann solver ( |rlarten et al. 198^ ; 
Einfeldt 1988). These fluxes are obtained independently for 
each direction and the time update of the state-vector U 
is done simultaneously using a method of lines in combina- 
tion with a second-order (in time) conservative Runge-Kutta 
scheme. Moreover, in order to set up a family of local Rie- 
mann problems at every cell-in terface we use a piecewise 
linear reconstruction procedure (van Leer 197£) which pro- 
vides second-order accuracy in space. 



4.2 Grid and boundary conditions 

We use a computational grid of 300 x 100 zones in the radial 
and angular direction, respectively. The grid is logarithmi- 
cally spaced in the radial direction. The innermost radius 
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is located at r m i n = 2.1. The location of the maximum 
radius r max depends on the particular model under study. 
For the stationary models presented in Section [s| we have 
r max = 35. Correspondingly, for the simulations of the run- 
away instability, the radial grid extends to a sufficiently large 
distance in order to ensure that the whole disc is included 
within the computational domain. The particular values of 
r max are displayed in Table ^ below. The typical width of 
the innermost cell, where we have the highest resolution, is 
Ar ~ 1.9 x 10~ 2 . 

In the angular direction we use a finer grid within the 
torus and a much coarser grid outside. The angular zones 
are distributed according to the following law: 



1 <j < 



"3 
7T 



1(0.1 + 2.8^) f + l<i<2M + i 
-0.2 + 1.2^ ^ + l<j<M + l 



(63) 



Although the flows we are simulating have equatorial 
plane symmetry we extend the computational domain in the 
angular direction from to 7T. This allows us to measure the 
ability of the code in keeping a symmetric evolution. 

For the second-order numerical scheme we use we need 
to impose boundary conditions in two additional zones at 
each end of the domain. The boundary conditions are ap- 
plied to p, v r ,ve and v$ and they are as follows: at the inner 
boundary r m i n all velocities are linearly extrapolated to the 
boundary zones from the innermost zones in the physical 
grid. The density is however assumed to have zero gradi- 
ent across the inner boundary. The rest of thermodynamical 
quantities are computed using the polytropic EoS. At the 
outer radial boundary r- max all variables keep the constant 
initial values given by the choice of the particular disc solu- 
tion. Reflection boundary conditions are used at both poles 
(6 — and 7r), i.e., all variables are symmetric, except for 
Vg which changes sign. 



4.3 Additional aspects 

Since we are considering adiabatic evolutions, we only solve 
for the first four equations of system (^). The internal en- 
ergy (proportional to the rest-mass density) is obtained alge- 
braically using a polytropic EoS, p = ftp 7 , i.e., e = ^-riP 7 ^ 1 . 
After the time update of the conserved quantities, the prim- 
itive variables are recomputed. As the relation between the 
two sets of variables is not in closed algebraic form, the prim- 
itive variables are computed using the following procedure: 
The evolved quantities D and Si being known, we eliminate 
p and h from the definition of Si given by Eq. ( |47[ ) to express 
the norm S 2 = (ph) 2 F 4 v 2 of Si as a function of the Lorentz 
factor r only: 



S 2 (T) = D 2 I 1 + 



7-1 



_D\ 7-1 

r 



(r 2 - 1) 



(64) 



We solve this equation by an iterative Newton-Raphson al- 
gorithm. Once the Lorentz factor F is found, the other primi- 
tive variables are easily derived using the relations p = D/F, 
h = l+ ^y^P 7-1 and v, = Si/ (phT 2 ) . 

Finally, it is worth pointing out that in order to evolve 
the "vacuum" zones which lie outside the disc using a hydro- 
dynamics code, we adopt the following simple yet effective 
procedure. Before constructing the initial torus we build up a 



background spherical accretion solution of a sufficiently low 
density so that its presence does not affect the dynamics of 
the disc. This stationary solution is given by the relativistic 
ext ension of the spherical Bondi accretion solution derived 
by Michel (1972). This solution depends on the location of 
the critical point r c and of the density at this point p c , to- 
gether with the adiabatic exponent and polytropic constant 
of the EoS, which we chose the same as inside the torus. In 
our approach we chose the values of r c and p c (which is com- 
puted from r c with the condition that the flow is regular at 
the critical point) in order to impose that the maximum den- 
sity in the background spherical solution is 5.0 x 10 -6 times 
the maximum density at the center of the disc. By doing 
this we have checked that the rest-mass present in our back- 
ground solution is always negligible compared to the mass 
of the disc and that the associated mass flux corresponding 
to this spherical accretion is also negligible compared to the 
mass flux from the disc. We note that in the outermost part 
the values of the background density can be as low as 10 -8 
times the maximum density at the center of the disc. 



5 SIMULATIONS OF STATIONARY MODELS 

As mentioned in the previous section our hydrodynamics 
code has been used previously in a number of relativis- 
tic wind accretion simulations. However, in order to test 
the ability of the code when dealing with accretion discs 
we have first considered time-dependent simulations of sta- 
tionary models. The aim of these simulations has been to 
find out whether the code is capable of keeping those mod- 
els in equilibrium during a sufficiently long period of time 
(much larger than the rotation period of the disc) . In order 
to do so we have considered the same stationary models that 
Igumenshchev & Beloborodov (1997) analyzed, in the limit 
of no black hole rotation. These four models are character- 
ized by I = 3.9136, a value in between the marginally sta- 
ble and marginally bound orbits, and an increasingly large 
value of the energy gap at the cusp, AWin = 0.02, 0.04, 
0.08 and 0.16. Similarly, the polytropic EoS has been chosen 
with an adiabatic index 7 = 4/3 and a polytropic constant 
k — 1.5 x 10 20 cgs. Furthermore, the mass of the black hole 
is kept constant throughout these test evolutions. 

As a representative example Figure ^ shows the mor- 
phology of the model with AW^n = 0.16. The cusp is located 
at r cusp = 4.197 and the center at reenter = 9.598. The grid 



extends to r max = 35. The top panel shows a gray-scale plot 
of the logarithm of the density for the initial model at t = 0, 
together with the corresponding velocity field. The bottom 
panel shows the same morphology at the final time t = 10 4 . 
This corresponds to about 53 dynamical timescales, choos- 
ing as a dynamical timescale the orbital period i or b = 2n/Q 
at r = reenter, which for this model is i or b = 184.9. The code 
was stopped after roughly 2 x 10 5 iterations with no signs of 
numerical instabilities present. One can clearly see in Fig- 
ure ^| that an accretion flow from the disc to the black hole 
appears in the inner region of the grid. This flow, however, 
becomes rapidly stationary (see below). In addition, for our 
particular choice of model parameters, the corresponding 
mass flux is very low. Therefore, except in the innermost 
region where the accretion flow develops, the morphology 
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Figure 3. Morphology of the inner part of the disc with param- 
eters I = 3.9136 and AlTi n = 0.16. Top: initial state; Bottom: 
final state at t ~ 53t or b . The arrows are proportional to the com- 
ponents (SV I yjg r r, Sg I y/gee^) °f the momentum and are plotted 
only in the region where p > 0.05 p„ lax . The black hole is repre- 
sented by the black circle and the exterior circle around it marks 
the location of the inner boundary of the grid. 



of the disc remains essentially unchanged during the whole 
evolution from t = to t = 53t or b- 

There are some additional noteworthy issues concern- 
ing this figure: first, it clearly shows the ability of the code 
to keep the equatorial plane symmetry of the torus, even 
though the angular domain extends from to n. Second, 
the smoothness of the initial disc distribution in the grid 
is maintained during the whole evolution even close to the 



black h ole horizon. Finally, contrary to previous work ( Haw- 



ley et cl. 1984b; [gumenshchev & Beloborodov 1997) there 



are no hints of vortices developing in the flow. In Hawley 



[gu- 



et al. (i 984b) such vortices are associated to small poloidal 
velociti es and are not as noticeabl e as in the results of 
mensh ofrev fc Beloborodov (199^ ). The latter claim, 1 
ever, that such vortex motions are likely to be related to the 
choice of initial conditions, developing from initial perturba- 
tions close to the cusp and propagating outwards undamped. 

A more quantitative proof of the ability of the code 
in keeping the stationarity of this solution is provided in 
Figure ^. This figure shows the evolution of the mass ac- 
cretion rate (normalized to the Eddington value), ra/mEdd, 
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Figure 4. Time evolution of the mass flux for a stationary model 
with I = 3.9136 and AW = 0.02, 0.04, 0.08 and 0.16. The time 
is given in units of the orbital period at the center (see text), 
which here is t or b — 187 r g /c. The mass flux is normalized with 
the Eddington limit computed with Mbh = 1 Mq. For a differ- 
ent choice of polytropic constant k and black hole mass Mbh> 
the mass flux scales as k~ 3 Mbh for 7 = 4/3. The initial value 
of the mass flux is fixed by the spherical mass accretion rate as- 
sociated with the Bondi flow imposed in the low density regions 
outside the torus. Notice that after a transition phase lasting for 
~ 0.1 i or b the mass accretion rate rapidly tends, asymptotically, 
to the stationary values. 

as a function of t/t OT h- The mass flux is computed at the 
innermost radial point as 



m = 2-k 



-gDv r d6, 



(65) 



where the volume element for the Schwarzschild metric is 
given by y/—g = r 2 s'm8. The Eddington mass flux riiEdd = 
~ 1.4 io 17 *£bu. g/ s is computed for M BH = 1 M . 
Rescaling of m for different polytropic constant k and black 
hole mass Mbh is given by 

{rh/rh E dd) 1 /ki\" 3 M B hi 



(m/rfiEdd) 



Mbh 



(66) 



where 7 = 4/3 was assumed. After a transient initial phase 
the mass accretion rate is seen to rapidly tend, asymptoti- 
cally, to a constant value. The offset observed during the ini- 
tial phase corresponds to the spherical accretion mass flux 
associated with the particular background solution we use 
outside the torus. 

Next, in Figure ^| we plot the mass flux as a function 
of the energy gap AWin for the four stationary models we 
have considered. The values selected for the mass flux in each 
model are the asymptotic ones, obtained after the simula- 
tions have been evolved up to a final time t = 10 4 , roughly 
54 orbital periods. This plot allows us to check if the code is 
able to reproduce the analytic dependence given by Eq. ( |32[ ) . 
For a 7 = 4/3 polytrope the expected slope is 4 (dashed 
line). As the figure shows our results are in good agreement 
with this analytic prediction as well as with the numer ical 
results obtained by Igumenshchev & Beloborodov (1997) for 
the same models. 
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AW in / e* 

Figure 5. Mass flux m as a function of the energy gap A Win 
between the inner edge of the disc and the cusp. The dots indicate 
the asymptotic values we get in our simulations (see Fig. ^). The 
dashed line shows the slope which is expected from the theoretical 
prediction given by Eq. ([j^), which for 7 = 4/3 is 4. The plot 
corresponds to simulations with Mbh = 1 Mq and k = 1.5 X 
10 20 cgs. 



6 SIMULATIONS OF THE RUNAWAY 
INSTABILITY 

6.1 The physical origin of the instability 

The physical mechanism le ading; to the runaway ins tabil- 
ity has been explain ed by Abramowicz ct al. (1983) and 
Nishida ct al. (1996). The mass transfer from the disc to 



the black hole starts once the disc fills its Roche lobe (see 
Fig. |^a). From this moment any small perturbation allows 
the gas to flow through the cusp located at the inner edge 
of the disc. As a result the mass of the black hole increases, 
the equipotential surfaces move, and the radial location of 
the cusp changes. The disc has to find a new equilibrium 
configuration: one possibility is that the disc overflows its 
Roche lobe as is depicted in Fig. |^b. In this case the mass 
transfer speeds up, which leads to the runaway instability. 
Alternatively, the disc may contract inside its Roche lobe 
(see Fig. ^c) . In this case the mass transfer slows down. The 
mass flux is self-regulated by this process and the accretion 
is stable. 

However, for a Schwarzschild black hole and a con- 
stant angular momentum distribution in the disc, most discs 
are unstable (Abramowicz et al. 198S) (see also Table [lj). 
This is illustrated in Fig. Q for the particular case where 
the black hole mass is A/bh = 2.44 Mq, the disc mass is 
Md = 0.36 Mq, the adiabatic index is 7 = 4/3 and the poly- 
tropic constant is k = 4.76 x 10 14 cgs (which corresponds to 
degenerate relativistic electrons with Y c = 0.5 electrons per 
nucleon) . In this figure the evolution of models with different 
mass is plotted against the mass transfered from the disc to 
the black hole. The initial disc is filling its Roche lobe. The 
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Figure 6. The physical origin of the runaway instability: 

the initial disc fills its Roche lobe (panel (a)). Mass transfer from 
the disc to the black hole occurs through the cusp, which is located 
at the Lagrange point L\. The mass of the black hole increases 
and the disc has to find a new equilibrium configuration with the 
new gravitational potential. There exist two possibilities: (i) the 
disc overflows its Roche lobe (panel (b)). This process speeds up 
the mass transfer and the disc becomes unstable, (ii) The disc 
contracts inside its Roche lobe (panel (c)), which forces the mass 
transfer to slow down, resulting on a stable disc. 
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Figure 7. The physical origin of the runaway instability: 

the evolution of models with different mass is plotted against 
the mass transfered from the disc to the black hole. The initial 
disc is filling its Roche lobe (indicated by the big dot symbol). 
Initially M BH = 2.44 M and Mn = 0.36 M Q . Upper panel: 
mass A/bh °f the black hole; Lower panel : mass Mo of the disc 
and mass Mj?J ax contained inside the Roche lobe (dashed line). 
As Mp > MJJ 1 ™, the disc overflows its Roche lobe which speeds 
up the mass transfer and leads to the runaway instability. 



location of both the initial disc mass and black hole mass is 
indicated by a big dot symbol in each panel. The upper panel 
shows the evolution of the mass of the black hole whereas 
the lower one shows the corresponding evolution of the disc 
and of M™, which is the maximum disc mass contained in- 
side the Roche lobe (dashed line). As Mn > Mg 1 ", the disc 
overflows its Roche lobe which speeds up the mass transfer 
and leads to the runaway instability. 



6.2 The sequence of stationary metrics 
approximat ion 

The increase of the mass of the black hole is the fundamen- 
tal process triggering the runaway instability. In order to 
properly take into account the dynamical evolution of the 
underlying gravitational field, one should solve the coupled 
system of Einstein and hydrodynamic equations. However, 
such a task has not yet been accomplished in the context 
of the runaway instability and, to some extent, it may be 
still far from the capabilities of current codes in numerical 
relativity. Numerical stability considerations, coming from 
both the coordinate singularity existing at the rotation axis 
(8 = 0, 7r), which spoils the long- term integration of the Ein- 
stein equations even in vacuum (Brandt & Seidel 1995), and 



from the mathematical formulation of the field equations 
themselves, make such a task very challenging. Numerical 
relativity codes evolving black hole spacetimes with perfect 
fluid matter are only becoming ava ilable very rec ently, both 
in axisymmetry ( Brandt et al. 200C ) as in full 3D ( Bhibata & 



Uryu 2000| ; |5hibata et al. 200u| ; [Alcubierre et al. 2000| ; |Fonl 
et al. 2001). Nevertheless, the resolution and the integration 



times required to study the time evolution of the runaway 
instability may be still too demanding for such relativistic 
codes which incorporate self-gravity. 

For all these reasons and for the complete lack of time- 
dependent simulations of the runaway instability in relativ- 
ity, we have adopted a simplified and pragmatic approach 
to the problem. In our procedure the spacetime metric is 
approximated at each time step by a stationary exact black 
hole metric of varying mass (and angular momentum in the 
case of a rotating black hole) . The mass M of the black hole 
necessary to compute the metric coefficients is increased at 
each time step At according to: 

M n+1 = M n + At ^ (g 7 ) 

where the mass flux at the inner radius of the grid is evalu- 
ated by the equation 



m n = 27rr, 



) sin 8jDi l 



■Lw- ( 68 ) 
As the mass of the black hole increases during the sim- 
ulations, the horizon moves outwards. To avoid the inner 
radius of the grid to become smaller than the radius of the 
growing horizon, we increase the index i m i n of the first radial 
zone when necessary, so that the condition rj min > 2M is al- 
ways respected. We notice that the black hole mass increases 
very slowly during the evolution which implies that the met- 
ric coefficients at any time step differ very little from the fi- 
nal values which would correspond to an exact Schwarzschild 
black hole of bigger mass but with no matter around. 

6.3 Initial state 

A given initial state of the black hole plus disc system is 
determined by five parameters: the mass of the black hole 
Mbh, the specific angular momentum in the disc I, the po- 
tential at the inner edge of the disc A Win, the polytropic 
constant k and the adiabatic index 7. The computing time 
needed for one hydrodynamical simulation is too large to al- 
low for a complete exploration of this parameter space. For 
this reason we focus on those models which are expected to 
be found in the central engine of GRBs. These systems are 
formed either after the coalescence of two compact objects 
or after the gravitational collapse of a massive star. 

6.3.1 Black hole mass and disc-to-hole mass ratio 

As it has been shown by numerical simulations using New- 
tonian and post-N ewtonian gravity, both the coalescence of 
two neutron stars ( Ruffert fc Janka 1999) and the mer ger of 
a black hole and a neutron star (Kluzniak & Lee 1998) lead 



to the formation of comparable systems, where the mass of 
the central black hole and the disc-to-hole mass ratio are 
respectively M BH ~ 2.5 M and M D /M BH ~ 0.04-0.08 for 
the first case, and Mbh ~ 3 M© and Mq/Mbh ~ 0.2 in the 
second case. More recently, the fully relativistic s imulations 
of binary ne utron star coalescence performed by Shibata & 
Uryu (200C) yield disc masses of ~ 0.05 — 0.1 A/» for corota- 



tional binaries and < 0.01M* for irrotational binaries, where 
M* is the total rest-mass of the system (typically ~ 2M©). 

The case of collapsars is more complex since, in par- 
ticular, the mass of the disc goes on increasing after the 



14 Font & Daigne 



Table 3. Initial models. The following parameters are listed: mass of the black hole Mbh> disc-to-hole mass ratio Mq/Mbh, specific 
angular momentum in the disc I, potential barrier at the inner edge AWi n , mass flux in the stationary regime m s t a t, minimum and 



maximum radii of the grid, r m i n and r max , radius of the cusp ?"cusp 5 radius of the center 7*0 



(all radii are in units of the gravitational 



radius r g ), and orbit al p eriod at the center of the disc t or b- The last column lists the timescale associated with the runaway instability 



as defined in Section 6.4. In all cases, the EoS parameters are k = 4.76 x 10 14 cgs and 7 = 4/3. 



Model 




M D /M BH 


I 




m s tat 




f*max 


f*cusp 


^"center 


^orb 


irun/^orb 




CM©) 








(M /s) 








(gco/ms) 






la 


2.5 


1. 


3.9325 


0.005 


0.090 


2.1 


85.* 


4.1492 


9.7930 


193. / 2.37 


97. / 95. + 


2a 


2.5 


1. 


3.9085 


0.01 


0.28 


2.1 


85.* 


4.2105 


9.5455 


185. / 2.28 


38. / 36. t 


3a 


2.5 


1. 


3.8564 


0.02 


2.1 


2.1 


85.* 


4.3639 


8.9919 


169. / 2.09 


10. / 9.2t 


4a 


2.5 


1. 


3.7255 


0.04 


34. 


2.1 


85.* 


5.0104 


7.3644 


126. / 1.55 


3.8 / 3.3 1 " 


lb 


2.5 


0.1 


3.8749 


0.005 


0.032 


2.1 


37. 


4.3058 


9.1915 


175. / 2.16 


140. 


2b 


2.5 


0.1 


3.8459 


0.01 


0.17 


2.1 


37. 


4.3990 


8.8768 


166. / 2.05 


64. 


3b 


2.5 


0.1 


3.7798 


0.02 


1.9 


2.1 


37. 


4.6698 


8.1077 


145. / 1.79 


12. 


lc 


2.5 


0.05 


3.8798 


0.001 


0.21 


2.1 


32. 


4.2911 


9.2438 


177. / 2.17 


110. 



* In models la to 4a the grid consists of a first grid from r = 2.1 to r = 28 and a second grid from r = 28 to r = 85 to avoid Ar at 
the inner radius becoming too large. 

' The second estimate of t run given for models la to 4a corres pon ds to simulations where the mass of the black hole starts to increase 
only once the stationary regime has been reached (see Section p. 4). 



formation of the black hole, as matter from the outer parts 
of the star is still infalling. T he si mulations perform ed by 
MacFadyen fc Woosley (1999|) and |Aloy ct al. (200(]|) start 



from tr e 14 Mq helium core of a rotating 35 M main se- 



quence star which collapses to produce a central black hole 
surrounded by a disc. When this systems enters a quasi- 
steady state, the mass of the black hole is 2-3 Mq and the 
disc-to-hole mass ratio is typically about 0.001-0.01. 

Taking into account these various results, we have de- 
cided to fix the mass of the black hole to Mbh = 2.5 Mq and 
to adjust the angular momentum I to get realistic disc-to- 
hole mass ratios. We have considered three possible ratios: 
0.05, 0.1 and 1. Such values are very close to what is ob- 
tained in the simulations of binary coalescence and not too 
far from the results of the simulations of collapsars. 



6.3.2 Equation of state 



We fix the adiabatic index to 7 
-.15 



4/3 and the polytropic 
constant to k = 1.2 x 10 la F c 4/3 with Y c = 0.5. This cor- 
responds to an EoS dominated by the contribution of rela- 
tivistic degenerate electrons (the typical density in the disc 
is ~ 10 n -10 12 g/cm 3 ). We note that such simplified EoS 
is nevertheless adequate to our purpose since the work of 
Nishida & Eriguchi (1996) showed that, for stationary mod- 
els, the effects of realistic EoS on the stability of constant 
angular momentum discs is negligible, the discs being un- 
stable in all cases. 



and 10 4 M /s. On the other hand, |Kluzniak fc Lee (199E 



find comparable values in the case of a neutron star - black 



hole merger. In t heir simulations of collapsars, MacFadyen 



fc Woosley (1999, ) find a typical mass flux of 0.6-0.8 M /s. 

Such mass fluxes are many orders of magnitude larger 
than the Eddington limit, which is 1.2 x 10~ 16 M /s for 
a Mbh = 2.5 M . However such very high mass fluxes are 
precisely what is required to explain the observed luminosity 
of GRBs, which is typically L 1 = 10 51 L51 erg/s in gamma- 
rays. If this radiation is due to internal shocks propagating 
within an ultra-relativistic wind, the kinetic energy flux of 
this wind is Lkm = L 7 // T , where / 7 is the efficiency of the 
kinetic energy to radiation conversion. If one assumes that 
the production of the relativistic wind is accretion-powered 
with an efficiency /ace, then the mass flux is given by 



1 L v 



°- 2 (mHc&r* 1 *®/' 



.A 



/a, 



(69) 



This estimate is of course no longer relevant if the main 
energy reservoir powering the burst is the rotational energy 
of the black h ole, which can be extracted by the Blandford- 
Znajek effect ([Blandford fc Znajek 19771). 



Taking into such high mass fluxes we have considered 
the following values of AW in : 0.001, 0.005, 0.01, 0.02 and 
0.04 so that the mass flux of our initial models in the sta- 
tionary regime spans ~ 0.03-30 M /s. From all the above 
considerations we have prepared eight initial models which 
are summarized in Table []. 



6.3.3 Mass flux 

Once Mbh, Md/Mbh, k and 7 have been specified, the only 
remaining parameter is the potential at the inner edge of the 
disc, AWin. We choose its value so that the corresponding 
mass flux in the stationary regime (as described in Section |)|) 
explores a realistic range. In the si mulations of binary neu - 
tron star coalescence carried out by Ruffert fc Janka (1999 ) , 
the mass accretion rate of the black hole varies between 1 



6.4 Results 

Each of our eight initial models is evolved twice: in the first 
series of runs we keep constant the mass of the black hole 
while in the second set of evolutions such mass is allowed 
to increase according to the law specified in the preceding 
section, Eqs. (^) and (jssj) . For the first type of runs the 
mass flux rapidly reaches a stationary regime (like in the 
models discussed in Section ^) which can be maintained for 
as long as desired. In practice the final time corresponds to 
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Figure 8. Time-evolution of the mass flux for models (la) to 
(4a) of Table |^. The solid lines correspond to evolutions in which 
the black hole mass incre ases with time according to the proce- 
dure explained in Section \>.2\ For comparison, the mass flux in 
the stationary regime (when the mass of the black hole is kept 
constant) is also plotted using dashed lines. As expected, for a 
black hole of growing mass the accretion process becomes rapidly 
unstable. Notice how the mass flux diverges. 



Figure 9. Time-evolution of the mass flux for models (lb) to 
(3b) and model (lc) of Table j?| As in Fig. |j, the solid lines corre- 
spond to evolutions in which the black hole mass inc rea ses with 
time according to the procedure explained in Section and the 
dashed lines correspond to evolutions in which the mass of the 
black hole is kept constant. As expected, for a black hole of grow- 
ing mass the accretion process becomes rapidly unstable. Notice 
how the mass flux diverges. 



many orbital periods of the disc. The mass flux at this stage 
has the value reported as rhstat in Table [| On the other 
hand, when the mass of the black hole increases (second 
series of runs), the time evolution of the system changes 
dramatically and the runaway instability appears. 

The dramatic differences between both series of simula- 
tions are depicted in Figures ^-[nj Figures |^ and show the 
time evolution of the mass accretion rate for models (la) to 
(4a) and (lb) to (3b) plus (lc), respectively. They include, 
both, the stationary and the unstable cases. The overall be- 
haviour found in both figures is similar despite the existing 
different mass ratio between the black hole and the disc, 1 in 
all curves of Fig. 8 and 0.1 and 0.05 in those of Fig. 9. The 
different mass ratio only affects the mass flux and (weakly) 
the time in which the instability appears. So, in models la- 
beled 'b' and 'c' the instability takes somewhat more time 
to appear than in the corresponding models 'a' with the 
same energy gap (see below). As a result, models labeled 'b' 
and 'c' need more computational time to be fully evolved. 
This fact sets important technical restrictions when trying 
to evolve models with smaller mass ratio for which neglect- 
ing the disc self-gravity would be more justified. Short-term 
evolutions of a model with a 0.01 disc-to-hole mass ratio 
show that the instability is indeed present but simply takes 
longer to grow. 

In Figures ^| and ^| we see that at early times the evo- 
lution of the mass flux for each pair of models is exactly 
the same, irrespective of the increase of the black hole mass 
being taken into account or not. However, whereas the mod- 
els with a constant black hole mass reach a quasi-stationary 



regime with a constant mass flux, the corresponding mass 
accretion rate for those models with an increasing black hole 
mass goes on increasing after having reached the 'stationary' 
value. Furthermore, the time derivative of the mass flux also 
increases, which implies that the mass flux diverges rapidly. 
For all unstable models computed the mass flux is already 
several orders of magnitude larger than the stationary value 
when the calculation is stopped. (4 orders of magnitude for 
Model (la)!, see Fig. ^). The divergence found in the mass 
accretion rate is a clear manifestation of the runaway in- 
stability at work, which leads ultimately to the complete 
destruction of the disc. 

In Figs. ^| and [H] we plot the time-evolution of the 
black hole mass and the disc mass for the same set of mod- 
els displayed in Figs. ^ and ^, respectively. The sudden loss 
of the mass of the disc at late times is reflected on the cor- 
responding rapid increase of the mass of the black hole. As 
an example, for Model (2a), at t ~ 40 t or b the black hole 
has almost doubled its mass (Mbh ~ 4.7 Mq) and, corre- 
spondingly, the mass of the disc has decreased from 2.5 Mq 
to roughly 0.3 Mq. Note that since models 'b' and 'c' have 
a much smaller disc mass than models 'a' the growth of 
the black hole mass in Fig. [ll] is not as clearly visible as in 
Fig. 0. 

The morphology changes that the unstable system un- 
dergoes are shown in Figure [l2] for a representative case 
(Model 3a). The evolution is qualitatively similar for all 
models. In this figure we show eight snapshots of the time- 
evolution from t = to t = 11.8 t or b- The variable plotted 
in the figure is the rest-mass density. The contour levels are 
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Figure 10. Time-evolution of the mass of the black hole and of 
the mass of the disc for models (la) to (4a) listed in Table |5| The 
sudden appearance of the instability is noticeable. The increase of 
the black hole mass is directly associated with the corresponding 
decrease in the mass of the disc. 
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Figure 11. Time-evolution of the mass of the black hole and 
of the mass of the disc for models (lb) to (3b) and model (lc) 
of Table ji| Since the mass of the disc is now much smaller than 
in models labeled 'a' the growth of the black hole mass is not as 
clearly visible as in Fig. llOl The sudden decrease of the mass of 
the disc is however noticeable. 



linearly spaced with Ap = 0.1 where p° is the maxi- 
mum value of the density at the center of the initial disc. In 
Fig. ^| one can clearly follow the transition from a quasi- 
stationary accretion regime (panels (1) to (5)) to the rapid 
development of the runaway instability (panels (6) to (8)). 
At t = 11.80 t or b, the disc has almost entirely disappeared 
inside the black hole whose size has noticeably grown. From 
the numerical point of view, and as already pointed out in 
Section |^ when describing stationary models, the flow solu- 
tion remains considerably smooth even though the evolution 
is now dynamic. The equatorial plane symmetry is main- 
tained during the whole evolution with no sign of numerical 
asymmetries as well as no vortices appearing inside the disc. 

Correspondingly, Figure [l3] shows the velocity field for 
model (3a) at t — 10.70t or b, associated with snapshot (7) in 
Fig. [l2]). This figure shows that the disc is falling radially 
on to the black hole with no signs of vortices and circulation 
patterns developing. 

An interesting information which our hydrodynamical 
simulations provide is the timescale of the instability t Iun . 
We estimate this timescale as the time it takes for half of the 
mass of the disc to fall into the hole. The values of i rU n ob- 
tained for our 8 models are given in the last column of Table 
^[ Such values span the interval ~ 3.8-140t or b, which corre- 
sponds to very small durations (~ 6-300 ms). To check the 
quality of our definition of t Iun we have performed the fol- 



lowing test: for the four models of series 'a' we have carried 
out additional simulations in which the mass of the black 
hole starts to increase only once the stationary regime has 
been reached at time t ~ to- The corresponding evolution of 
the mass flux in case (4a), for which to — 7.9 t OI b, is plotted 
in Fig. |l4| To compare more easily the timescale associated 
with the runaway instability in the two cases, we have plot- 
ted in Fig. [Hi] the evolution of the mass flux as a function of 
t — to ■ Again we see that in all cases the runaway instability 
appears immediately resulting in the rapid disappearance of 
the disc. However the precise comparison between our usual 
series of runs (where to = 0) and the modified ones leads to 
the conclusion that the timescale of the runaway instability 
is a bit overestimated in the first case. The new values of 
the timescale for series 'a' are reported in the last column 
of Table |. 

In Fig. [l(j we plot the timescale t Iun as a function of the 
mass flux in the stationary regime m sta t- To be able to de- 
rive an empirical law for the timescale of the instability one 
should consider a much larger sample of models. Neverthe- 
less, two clear tendencies can already be extracted from this 
figure: (i) the timescale depends weakly on the disc-to-hole 
mass ratio; (ii) the runaway instability occurs faster when 
the initial mass flux (stationary value) is larger, following 
approximatively t Tun oc m~ a . With the value of a = 0.9 ob- 
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* / rg 

Figure 12. Time evolution of the unstable model (3a): contour levels of the rest-mass density p plotted at irregular times from t = to 
t = 11.80 t or b once the disc has almost been entirely destroyed. The density levels are linearly spaced with Ap = 0.1 p®. The density p® 
takes the value at the center of the initial disc (marked by a dot in panel (1)). The most exterior contour corresponds to p = 0.1 p^ax- 
The entire disc is filled in grey color. From panels (1) to (3) the disc is very close to a stationary regime and it is almost not evolving. 
The runaway instability develops from panels (4) to (8), most noticeable in the last three panels from time t — 9.47 t rh 

to t = 11.80 i orh . 

The increase in the central density and the infall of the disc to the black hole are well visible, as well as the associated growth of the 
black hole. 
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Figure 13. Unstable Model (3a): velocity field at t = 10.70 f orb 
(corresponding to snapshot (7) in Fig. |l2|). The arrows are propor- 
tional to v' and are plotted only in the region where p > 0.1 p^ax- 
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Figure 14. Time-evolution of the mass flux for model (4a) of 
Table jj| As in Figure ^, the dashed line corresponds to the case 
where the mass of the black hole is kept constant whereas the 
thin solid line corresponds to the case where the mass of the 
black hole increases (from t = 0). In addition, the thick solid 
line corresponds to a third case in which the mass of the black 
hole starts to increase only at time to = 7.9 t or b (indicated by a 
vertical arrow), i.e. once the mass flux has reached its 'stationary 
value'. In this case the instability shows up immediately. However, 
the timescale is the same than in the previous run (thin solid line; 
notice the logarithmic scale for t). 



tained for series 'a' (where we have used the result of the 
modified runs to get a more accurate estimate of t run ), we 
can infer that, for all cases, the disc is destroyed in a dura- 
tion never exceeding 1 s for a large range of accretion mass 
fluxes, rhstat > 10~ 3 M /s. 




(t-y/t„ b 

Figure 15. Time-evolution of the mass flux as a function of t — to 
for models (la) to (4a) of Table ^, where to is the time at which 
the mass of the black hole starts to increase. As in Figure ^, the 
dashed lines correspond to the case where the mass of the black 
hole is kept constant whereas the thin solid lines correspond to 
the case where the mass of the black hole increases (to = 0). The 
thick solid lines show the evolution when the mass of the black 
hole starts to increase only from time to, equal to 1000, 200, 24 
and 7.9 t or b for models (la) to (4a) respectively, i.e. once the mass 
flux has reached its 'stationary value'. In such case, the runaway 
instability appears earlier, the effect being more important for 
models (la) and (2a). 



7 CONCLUSIONS 

We have presented results from a numerical study of the 
runaway instability of thick discs around black holes. In 
this study we have carried out a comprehensive set of time- 
dependent simulations aimed at exploring the appearance 
of the instability. In order to do so we have used a fully 
relativistic, axisymmetric hydrodynamics code. The general 
relativistic hydrodynamic equations have been formulated as 
a first-order, flux-conservative hyperbolic system and solved 
using a suitable Godunov-type scheme. Among the simplify- 
ing conditions we have assumed a constant angular momen- 
tum disc around a Schwarzschild (nonrotating) black hole. 
The self-gravity of the disc has been neglected and the evo- 
lution of the central black hole has been assumed to be that 
of a sequence of exact Schwarzschild black holes of varying 
mass. 

We have found that by allowing the mass of the black 
hole to grow the runaway instability appears on a dynamical 
timescale. The mass flux diverges and the disc entirely falls 
into the hole in a few orbital periods (1 — > 100). Therefore, 
the appearance of the runaway instability in constant angu- 
lar momentum discs found in our simulations is in complete 
agreement with previous estimates from stat ionary models 



(Abramowicz et al. 1983; Nishida et al. 1996). Our simula- 



tions provide the first estimation of the timescale associated 
with the runaway instability. For a black hole of 2.5 Mq and 
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gular momentum in the disc (i.e. increasing outwards) has 
also not been considered yet. 

The first two points in the above list cannot be im- 
proved without solving the coupled system of Einstein and 
hydrodynamic equations on black hole spacetimes. Although 
important advances in the field of numerical relativity this 
task is still challenging. The other items, however, can be 
more easily improved and work in this direction will be pre- 
sented in subsequent investigations. In particular, we will 
present in a forthcoming paper the effect of a non-constant 
angular momentum in the disc. Such a distribution of angu- 
lar momentum is believed to suppress the runaway instabil- 
ity a ccording to previous studies in a stationary framework 
(e.g. Daigne fc Mochkovitch (1997 )). This last point - and 
the very existence of the runaway instability itself - is very 
important in the context of the most discussed scenario for 
GRBs. In the standard model the central engine responsi- 
ble for the highly energetic emission is a thick disc orbiting 
a stellar mass black hole, with a high accretion mass flux. 
The lifetime of this system must necessarily be larger than a 
few seconds to explain the observed durations of the bursts. 
Our results show that it would be absolutely excluded if the 
runaway instability occurs. 



Figure 16. The time scale of the runaway instability t run (as 
defined in Section 6.4 ) is plotted as a function of the mass flux 



"istat in the 'stationary regime' for each of the models listed in 
Table |. For models (la) to (4a) (M D /M B H = !•), we use the 
more accurate estimate of t Tun obtained when the mass of the 
black hole starts to increase once the stationary regime has been 
reached. The dashed line corresponds to the best fit 

^run ^otat 

for this series. We find a = 0.9. 



disc-to-hole mass ratios between 1 and 0.05 this timescale 
never exceeds 1 s for a large range of mass fluxes and it is 
typically about 50 ms. We have found that the dependence 
of the timescale on the disc-to-hole mass ratio is weak and 
that the runaway instability occurs faster the larger it is the 
initial mass flux (stationary regime) from the disc to the 
black hole. 

We note that our study has been restricted to a poly- 
tropic gas, with a particular choice of n and 7 corresponding 
to a gas of degenerate relativistic electrons. We are aware of 
th e over-simplification of su ch an EoS. However, the work 
of Nishida fc Eriguchi (1996 ) has shown that the conclusion 
of Nishida et al. (1991:) (where stationary models were built 
in a fully relativistic computation including the self-gravity 
of the disc) is not modified when using a realistic EoS: con- 
stant angular momentum discs are unstable. Therefore, we 
believe that our results would not be strongly modified if we 
were using a more elaborate description of the matter. 

To close our investigation we notice that there are four 
important limitations in our study: (i) it is difficult to check 
the validity of our simple-minded approach to incorporate 
the effect of the black hole mass increase; (ii) the self-gravity 
of the disc has not been included. Studies based on sequences 
of equilibri um configurations have shown that it favo rs the 
instability flNishida et al. 1996); [Masuda et al. 1998[); (iii) 



the rotation of the black hole and the possible increase of its 
spin due to the transfer of angular momentum (associated 
with the transfer of mass) is not yet included in our current 
model; (iv) the case of a more realistic distribution of an- 
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